function [U,bf,p,X] = fsb_conv_onsets(inparg)

% FSB - DEV: convolve time-jittered files for sandbox
% with SPM functions to hrf, optionally with time and dispersion
% derivatives.
%
% EXAMPLE:
% [U,bf,p,X, hemodynamics] = fsb_conv_onsets(inparg)
%
% INPUT:
% inparg: struct containing the fields
% evt_onsets : event onsets e.g. (33.5 46.5 48) and evt_dur
% evt_names : event names e.g. 'Eyes', 'Faces', 'Houses';
% evt_dur : event durations e.g. (3 3 2);
% TR : TR;
% volumes : number of volumes in scan;
% deriv : set this to one if you want to have the derivatives, too
%
% OUTPUT:
% U
% bf
% p
% X: hemodynamics array for Sandbox
%
% CALLED BY:
% FSB.m
% fsb_get_hemodynamics_er.m
% fsb_get_hemodynamics_er_baseline.m
% fsb_get_hemodynamics_blocked.m
% fsb_hrf_auto_adjust.m
%
% NOTES:
% Code adapted from SPM2 functions 
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% first part: spm2/spm_get_bf.m
% second part: spm2/spm_get_ons.m
% 
% this function relies on the following SPM2 functions: 
% spm_hrf.m and spm_Gdpf.m for hrf calculation
% spm_orth.m
%
% SPM is free but copyright software, distributed under the terms of the GNU
% General Public Licence as published by the Free Software Foundation
% (either version 2, as given in file spm_LICENCE.man, or at your option,
%
% this function assumes that no slice correction is done and shifts the hemodynamics
% predictor accordingly
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $ Revision 1.0
%
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

T     = 16;
fMRI_T  = 16;
fMRI_T0 = 10;

dt=(1/T)*inparg.TR;
evt_onsets = inparg.evt_onsets;
%evt_names = inparg.evt_names;
evt_dur = inparg.evt_dur;
nscans = inparg.volumes;

noonset =[];

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check input format for evt_dur
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

for i=1:length(evt_onsets);
    %
    try
        evt_durs{i} = evt_dur{i}; % if evt_dur is a cell
    catch
        evt_durs{i} = evt_dur(i)*ones(1, length(evt_onsets{i})); % if evt_dur is a vector
    end
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % insert a shift here to account for lack of slice time correction
    % (onset half volume too late)
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    evt_onsets{i} = evt_onsets{i}-0.5;
end;

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Set up basis functions with canonical hemodynamic response function
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[bf p]         = spm_hrf(dt);

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate derivatives if so desired
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if inparg.deriv ==1;
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % time derivative
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    dp     = 1;
    p(6)   = p(6) + dp;
    D      = (bf(:,1) - spm_hrf(dt,p))/dp;
    bf     = [bf D(:)];
    p(6)   = p(6) - dp;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % dispersion derivative
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    dp    = 0.01;
    p(3)  = p(3) + dp;
    D     = (bf(:,1) - spm_hrf(dt,p))/dp;
    bf    = [bf D(:)];

end

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Orthogonalize and fill in basis function structure
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bf  =  spm_orth(bf);

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create model
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

for i = 1:length(evt_onsets);
   % for i = 1:size(evt_onsets,2);

    if ~isempty(evt_onsets{i}) % Check if onsets exist for condition
        noonset(i) = 0;

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % onsets
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ons = evt_onsets{i};
        ons = ons(:);
        u     = ons.^0;

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % durations
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        dur = evt_durs{i};
        dur = dur(:);

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % peri-stimulus times {seconds}
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        pst   = [1:nscans]*T*dt - ons(1)*inparg.TR;
        for j = 1:length(ons)
            w      = [1:nscans]*T*dt - ons(j)*inparg.TR;
            v      = find(w >= -1);
            pst(v) = w(v);
        end
        
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Create dummy information
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        xP.name  = 'none';
        xP.h     = 0;

        Uname     = inparg.evt_names{1,i};
        
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~        
        % interaction with causes (u) - 1st = main effects
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        u     = ons.^0;
        for q = 1:length(xP)
            xP(q).i = [1, ([1:xP(q).h] + size(u,2))];
            for   j = 1:xP(q).h
                u   = [u xP(q).P.^j];
                str = sprintf('%sx%s^%d',Uname{1},xP(q).name,j);
                Uname{end + 1} = str;
            end
        end

        
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % orthogonalize inputs
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        u          = spm_orth(u);

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % and scale so sum(u*dt) = number of events, if event-related
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if ~any(dur)
            u  = u/dt;
        end
        
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % create stimulus functions (32 bin offset)
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ton       = round(ons*inparg.TR/dt) + 32;			% onsets
        tof       = round(dur*inparg.TR/dt) + ton + 1;			% offset
        sfmat        = sparse((nscans*T + 128),size(u,2));

        for j = 1:length(ton)
            sfmat(ton(j),:) = sfmat(ton(j),:) + u(j,:);
            sfmat(tof(j),:) = sfmat(tof(j),:) - u(j,:);
        end

        sfmat        = cumsum(sfmat);					% integrate
        sfmat        = sfmat(1:(nscans*T + 32),:);				% stimulus

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % place in ouputs structure
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        U(i).u    = sfmat;			% - stimulus function matrix

    else

        noonset(i) = 1;

    end
end

if length(noonset)>0; % if there is at least one onset vector
    xx = zeros(length(sfmat),1);
    for x = 1: length(noonset);
        if noonset(x) ==1;
            U(1,x).u = xx;
        end
    end

    X     = [];

    for i = 1:length(U)

        for nscans = 1:size(U(i).u,2)
            for p = 1:size(bf,2)

                x      = U(i).u(:,nscans);
                d      = 1:length(x);
                x      = conv(full(x),bf(:,p));
                x      = x(d);
                X      = [X x];

            end
        end

    end

    X = X([0:(inparg.volumes - 1)]*fMRI_T + fMRI_T0 + 32,:);
end
end


